Author: Amanda Everitt
Began: 8/22/2018
Finished:: 8/23/2018
[Goal]
- Identify which cell types are present. We would expect most of these are neuronal.
[Results]
- 14933 genes after filtering
- 17396 cells after filtering
- Data was normalized (~log(cpm)), scaled, centered, and regressed (by nUMI, nGene, %mito)
- 12 clusters identified, 7 cell types identified
| Neuronal |
70.52% |
0,1,2 |
Nrgn, Cnih2 |
| Interneuron |
11.54% |
5,6,8 |
Dlx1, Dlx2, Arx, Gad2 |
| Astrocyte |
4.98% |
4 |
Aldh1l1, Aqp4, Gja1, Mfge8 |
| Microglia |
5.31% |
3 |
Aif1, C1qa, C1qb, Ctss |
| Endothelial |
3.55% |
7 |
Cldn5, Ctla2a, Igfbp7, Kdr |
| Pericyte |
2.14% |
10,11 |
Igf2, Col1a1, Col1a2, Dcn |
| Oligodendrocyte |
1.92% |
9 |
Gpr17, Olig2, Pdgfra, Itpr2 |
1. Load Data
load(paste0(wd, "/outputs/output_01/experiment.aggregate.filtered.Rdata"))
#Normalize -> Log(CPM)
experiment.aggregate <- NormalizeData(object= experiment.aggregate, normalization.method = "LogNormalize", scale.factor = 10000)
2. Exploratory look at data
genes_to_plot = c("Tbr1", "Fezf2","Bcl11b","Bcl11a","Grin2b","Foxp1","Nrgn","Apoe")
plots_list <- list()
for (i in genes_to_plot){
a <- as.data.frame(Matrix::colSums(experiment.aggregate@data[i, , drop = FALSE]))
a$sample_type <- "WT"
a[grep("NULL", rownames(a)),]$sample_type <- "NULL"
a[grep("HET", rownames(a)),]$sample_type <- "HET"
colnames(a) <- c("count","sample_type")
p1<- ggplot(a, aes(x=count, fill=sample_type)) +
geom_density(alpha= 0.3) +
scale_x_continuous(name="Normalized counts per cell") +
ylim(0,1) + ylab("Density") +
ggtitle(paste0("Density of ", i," expression"))
plots_list[[i]] <- p1
}
marrangeGrob(grobs=plots_list, nrow=2, ncol=2)


GenePlot( experiment.aggregate, "nUMI", "nGene", cex.use = 0.5)

3. Identify variable genes
experiment.aggregate <- FindVariableGenes(
object = experiment.aggregate,
mean.function = ExpMean,
dispersion.function = LogVMR,
x.low.cutoff = 0.125,
x.high.cutoff = 4,
y.cutoff = 0.5, do.plot=F)
length(experiment.aggregate@var.genes)
## [1] 528
4. Linear Dimensional Reduction (sequencing depth, and % mito)
#Scale and Regress data on sequencing depth and % mito
experiment.aggregate <- ScaleData(
object = experiment.aggregate, #will populate new slot "scale data", does not change "data" slot
do.scale = TRUE,
do.center = TRUE,
vars.to.regress = c("nUMI", "nGene","percent.mito"))
## Regressing out: nUMI, nGene, percent.mito
##
## Time Elapsed: 1.24776516755422 mins
## Scaling data matrix
#Build PCs
experiment.aggregate <- RunPCA(
object = experiment.aggregate,
pc.genes = experiment.aggregate@var.genes, #use variable genes
pcs.compute = 20,
maxit = 500)
## [1] "PC1"
## [1] "Tubb2b" "Stmn1" "Stmn2" "Tubb3" "Rtn1" "Fabp7" "Gap43"
## [8] "Map1b" "Calm2" "Ncam1" "Nnat" "Sox11" "Snrpn" "Nrgn"
## [15] "Gm3764" "Meis2" "Mt3" "Igfbpl1" "Gria2" "Nrxn3" "Tubb5"
## [22] "Sox2" "Ccnd2" "Ckb" "Dbi" "Atpif1" "Dpysl2" "Ptprz1"
## [29] "Dlx1" "Dlx2"
## [1] ""
## [1] "Cldn5" "Esam" "Igfbp7" "Flt1" "Vwa1" "Pglyrp1"
## [7] "Itm2a" "Slc38a5" "BC028528" "Egfl7" "Ctla2a" "AU021092"
## [13] "Slc22a8" "Ecscr" "Slc16a1" "Cd34" "Kdr" "Eng"
## [19] "Pltp" "Cd93" "Fn1" "Ctsh" "Lmo2" "Slc40a1"
## [25] "Gng11" "Ramp2" "Car4" "Slc7a1" "Klf2" "Slc7a5"
## [1] ""
## [1] ""
## [1] "PC2"
## [1] "Cldn5" "Esam" "Vwa1" "Slc38a5" "Igfbp7" "Flt1"
## [7] "Pglyrp1" "Ctla2a" "Itm2a" "AU021092" "Slc22a8" "Kdr"
## [13] "Slc16a1" "Mfsd2a" "Slc7a1" "Car4" "Spock2" "Fn1"
## [19] "Ly6c1" "Ramp2" "Eng" "Lsr" "Cd93" "Slc7a5"
## [25] "Gng11" "Cd34" "Slco1c1" "Tfrc" "Ets1" "Htra3"
## [1] ""
## [1] "C1qc" "C1qa" "C1qb" "Csf1r" "Fcrls" "Fcer1g" "Ctss"
## [8] "Tyrobp" "Pld4" "Laptm5" "Fcgr3" "Trem2" "Cd68" "Cx3cr1"
## [15] "Ly86" "Ltc4s" "Lgmn" "Unc93b1" "P2ry12" "Hexb" "Lyz2"
## [22] "Spi1" "Gpr34" "Grn" "Aif1" "Ctsz" "Gmfg" "Ctsd"
## [29] "Hexa" "Olfml3"
## [1] ""
## [1] ""
## [1] "PC3"
## [1] "Stmn1" "Igfbpl1" "Tubb5" "Dlx1"
## [5] "Dlx6os1" "Dlx2" "Tubb3" "Stmn2"
## [9] "Sp8" "Calm2" "Ptma" "Atpif1"
## [13] "Sp9" "Nrxn3" "Dlx5" "Meg3"
## [17] "Cdca7" "Sox11" "Hmgn2" "Hist1h2ap"
## [21] "2810417H13Rik" "Map1b" "Gad2" "Tiam2"
## [25] "Hmgb1" "Arx" "Top2a" "Rnd3"
## [29] "Sox4" "Actb"
## [1] ""
## [1] "Apoe" "Slc1a3" "Ptprz1" "Plpp3" "Mfge8" "Ptn" "Bcan"
## [8] "Gja1" "Aqp4" "Pla2g7" "Atp1a2" "Ednrb" "Htra1" "Slc4a4"
## [15] "Sparcl1" "Fabp7" "Igfbp2" "Slc7a10" "Atp1b2" "Cd9" "Cspg5"
## [22] "Tril" "Gpr37l1" "Lfng" "Btbd17" "Cd63" "Clu" "Slc6a11"
## [29] "Dbi" "Vcam1"
## [1] ""
## [1] ""
## [1] "PC4"
## [1] "Gap43" "Nrgn" "Snrpn" "Calm1" "Calm2" "Ndufa4" "Uqcrq"
## [8] "Stmn2" "Usmg5" "Snhg11" "Pcp4" "Atpif1" "Dbi" "Ache"
## [15] "Cplx2" "Reln" "Opcml" "Syp" "Mt3" "Bhlhe22" "Fabp7"
## [22] "Tubb2b" "Spock1" "Sst" "Lrfn5" "Clstn3" "Scg2" "Tubb3"
## [29] "Vgf" "Stmn1"
## [1] ""
## [1] "Igfbpl1" "Top2a" "Sp8" "Zbtb20"
## [5] "Dlx2" "Nfib" "Hmgb2" "Dlx6os1"
## [9] "Cdca7" "Hist1h2ap" "Dlx1" "Xist"
## [13] "Nusap1" "Sp9" "Meis2" "2810417H13Rik"
## [17] "Ube2c" "Ccnb1" "Hmgn2" "Sox11"
## [21] "Tpx2" "Ccnd2" "Cenpf" "Cks2"
## [25] "Cdca8" "Birc5" "Cdca3" "H2afv"
## [29] "Cenpa" "Cdc20"
## [1] ""
## [1] ""
## [1] "PC5"
## [1] "Vtn" "Igf2" "Col1a2" "Cald1" "Col1a1" "Myl9"
## [7] "Ndufa4l2" "Dcn" "Emp3" "Rgs5" "Higd1b" "Nupr1"
## [13] "Rarres2" "Phlda1" "Fstl1" "Bgn" "S100a11" "Col4a1"
## [19] "Nbl1" "Nr2f2" "Ifitm3" "Serpinf1" "Zic1" "Cthrc1"
## [25] "Lgals1" "Eva1b" "Col4a2" "Cxcl12" "Ifitm2" "Gpx8"
## [1] ""
## [1] "Cldn5" "Slco1c1" "Slc38a5" "AU021092" "Ctla2a" "Slc7a5"
## [7] "Pglyrp1" "Slc16a1" "Mfsd2a" "Kdr" "Slc7a1" "Spock2"
## [13] "Vwa1" "Tfrc" "Ly6c1" "Lmo2" "Cd34" "Slc40a1"
## [19] "Cd93" "Egfl7" "Slc2a1" "Flt1" "Nampt" "Car4"
## [25] "Htra3" "Ptprz1" "Fli1" "Slc1a4" "Lsr" "Slc22a8"
## [1] ""
## [1] ""
PCAPlot(object = experiment.aggregate,dim.1 = 1, dim.2 = 2 )

VizPCA(object = experiment.aggregate, pcs.use=1:2)

PCHeatmap(
object = experiment.aggregate,
pc.use = 1:12,
cells.use = 500,
do.balanced = TRUE,
label.columns = FALSE,
use.full = FALSE
)

5. Determine statistically significant PCs
experiment.aggregate <- JackStraw(
object = experiment.aggregate,
num.replicate = 100,
num.pc = 20)
#Too slow above so loading here instead of regenerating
load(file=paste0(wd, "/", out_dir, "/experiment.aggregate.norm.RData"))
PCElbowPlot(experiment.aggregate,num.pc = 40)

JackStrawPlot(object = experiment.aggregate, PCs = 1:40, nCol = 5)
## Warning: Removed 16527 rows containing missing values (geom_point).

## An object of class seurat in project Siavash scRNA
## 14933 genes across 17396 samples.
percentVar <- experiment.aggregate@dr$pca@sdev^2/sum(experiment.aggregate@dr$pca@sdev^2)*100
d<- as.data.frame(percentVar)
d$PC <- as.factor(paste0("PC",rownames(d)))
positions <- d$PC
p<-ggplot(data=d, aes(x=PC, y=percentVar)) +
geom_bar(stat="identity",fill="steelblue") +
theme(axis.text.x=element_text(angle=90, hjust=1), legend.position="none") +
#geom_text(aes(label=round(percentVar,digits = 1)), vjust=1.6, color="white", size=3.5)+
xlab("Principle Component") +
ylab("Variance Explained (%)") +
scale_x_discrete(limits = positions)+
scale_y_continuous(breaks = seq(0, 20, by = 1))+
geom_hline(yintercept=2, col="red")
p

use.pcs = c(1:11) #Chose all PCs accounting for more than 2% of variance
experiment.aggregate <- FindClusters(
object = experiment.aggregate,
reduction.type = "pca",
dims.use = use.pcs,
resolution = seq(0.3,1.2,0.3),
print.output = FALSE,
force.recalc = TRUE,
save.SNN = TRUE
)
save("experiment.aggregate", file=paste0(wd, "/", out_dir, "/experiment.aggregate.norm.RData"))
use.pcs = c(1:11)
sapply(grep("^res",colnames(experiment.aggregate@meta.data),value = TRUE),
function(x) length(unique(experiment.aggregate@meta.data[,x])))
## res.0.3 res.0.6 res.0.9 res.1.2
## 12 15 18 20
6. Choose which resolution is best
experiment.aggregate <- SetAllIdent(experiment.aggregate, id = "res.0.3") #temporarily chose one
experiment.aggregate <- RunTSNE( object = experiment.aggregate, reduction.use = "pca", dims.use = use.pcs, do.fast = TRUE)
load(file=paste0(out_dir, "/experiment_aggregate.0.3.RData"))
TSNEPlot(object = experiment.aggregate.0.3, group.by="res.0.3", pt.size=0.5, do.label = T)

TSNEPlot(object = experiment.aggregate.0.3, group.by="res.0.6", pt.size=0.5, do.label = T)

#I chose resolution 0.3 here
experiment.aggregate.0.3 <- experiment.aggregate
table(experiment.aggregate.0.3@ident,experiment.aggregate.0.3@meta.data$orig.ident)
experiment.aggregate.0.3 <- BuildClusterTree(experiment.aggregate.0.3,
pcs.use = use.pcs, do.reorder = F,
reorder.numeric = F, do.plot=F)
save(experiment.aggregate.0.3, file=paste0(wd, "/", out_dir, "/experiment_aggregate.0.3.RData"))
TSNEPlot(object = experiment.aggregate.0.3, pt.size=0.5, do.label=TRUE)

TSNEPlot(object = experiment.aggregate.0.3, group.by="orig.ident", pt.size=0.5)

PlotClusterTree(experiment.aggregate.0.3)

7. Useful plots
#Silhoutte Lengths-- how close is this clustering?
pcas <- experiment.aggregate.0.3@dr$pca@cell.embeddings[, use.pcs]
my.dist <- dist(pcas)
my.clusters <- as.numeric(as.factor(experiment.aggregate.0.3@meta.data$res.0.3))
clust.col <- hue_pal()(13)
sil <- silhouette(my.clusters, dist = my.dist)
sil.cols <- clust.col[ifelse(sil[,3] > 0, sil[,1], sil[,2])]
sil.cols <- sil.cols[order(-sil[,1], sil[,3])]
plot(sil, main = paste(length(unique(my.clusters)), "clusters"),
border=sil.cols, col=sil.cols, do.col.sort=FALSE)

- Barplot of silhouette widths for cells in each cluster. Each cluster is assigned a colour and cells with positive widths are coloured according to the colour of its assigned cluster. Any cell with a negative width is coloured according to the colour of the cluster that it is closest to. The average width for all cells in each cluster is shown, along with the average width for all cells in the data set.
See if the possible outliter cells from scater should be removed (I decide not to)
load(paste0(wd, "/outputs/output_01/possible_outliers.Rdata"))
experiment.aggregate.0.3@meta.data$scater_outlier <- "NO"
experiment.aggregate.0.3@meta.data[rownames(experiment.aggregate.0.3@meta.data) %in%
possible_outliers,]$scater_outlier <- "YES"
table(experiment.aggregate.0.3@meta.data$scater_outlier)
##
## NO YES
## 16702 694
TSNEPlot(object = experiment.aggregate.0.3, group.by="scater_outlier", pt.size=0.5, do.return = TRUE)

8. Output Marker Genes for each cluster
markers_all.0.3 <- FindAllMarkers(object = experiment.aggregate.0.3, min.pct = 0.25,
thresh.use = 0.25, only.pos = TRUE,
test.use = "wilcox")
write.table(markers_all.0.3, file = paste0(wd, "/", out_dir, "/01_marker_genes.csv"), sep=",", col.names = NA)
save("markers_all.0.3", file=paste0(wd, "/", out_dir, "/markers_all.0.3.RData"))
load(paste0(wd, "/", out_dir, "/markers_all.0.3.RData"))
| Neuronal |
70.52% |
0,1,2 |
Nrgn, Rorb, Cnih2 |
| Interneuron |
11.54% |
5,6,8 |
Dlx1, Dlx2, Arx, Gad2 |
| Astrocyte |
4.98% |
4 |
Aldh1l1, Aqp4, Gja1, Mfge8 |
| Microglia |
5.31% |
3 |
Aif1, C1qa, C1qb, Ctss |
| Endothelial |
3.55% |
7 |
Cldn5, Ctla2a, Igfbp7, Kdr |
| Pericyte |
2.14% |
10,11 |
Igf2, Col1a1, Col1a2, Dcn |
| Oligodendrocyte |
1.92% |
9 |
Gpr17, Olig2, Pdgfra, Itpr2 |
my.col = c("lightblue", "red")
#Neuronal
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Nrgn","Rorb"), cols.use = my.col, nCol = 2, no.legend = F)

#Interneron
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Dlx1","Dlx2","Arx","Gad2"), cols.use = my.col, nCol = 2, no.legend = F)

#Astrocyte
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Aldh1l1","Aqp4","Gja1","Mfge8"), cols.use = my.col, nCol = 2, no.legend = F)

#Microglia
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Aif1","C1qa","C1qb","Ctss"), cols.use = my.col, nCol = 2, no.legend = F)

#Endothelial
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Cldn5","Ctla2a","Igfbp7","Kdr"), cols.use = my.col, nCol = 2, no.legend = F)

#Pericyte
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Igf2","Col1a1","Col1a2","Dcn"), cols.use = my.col, nCol = 2, no.legend = F)

#Oligo
FeaturePlot(experiment.aggregate.0.3, features.plot = c("Gpr17","Olig2","Pdgfra","Itpr2"), cols.use = my.col, nCol = 2, no.legend = F)

9. Relabel clusters with intial DEX results
#TSNE
current.cluster.ids <- c(0:11)
new.cluster.ids <- c("Neuronal", "Neuronal", "Neuronal","Microglia","Astrocyte",
"Interneuron","Interneuron","Endothelial","Interneuron",
"Oligodendrocytes","Pericyte","Pericyte")
experiment.aggregate.0.3@ident <- plyr::mapvalues(x = experiment.aggregate.0.3@ident,
from = current.cluster.ids, to = new.cluster.ids)
p1 <- TSNEPlot(object = experiment.aggregate.0.3, pt.size = 0.1, do.return=TRUE)
p1 <- p1 + ggtitle("tSNE by cell type")
p2 <- TSNEPlot(object = experiment.aggregate.0.3, group.by="orig.ident", pt.size = 0.1, do.return=TRUE, colors.use = c("#F8766D","#7CAE00","skyblue"))
p2 <- p2 + ggtitle("tSNE by genotype")
grid.arrange(p1,p2,ncol=2)

#Cluster Tree
suppressPackageStartupMessages(library(ape))
data.tree <- experiment.aggregate.0.3@cluster.tree[[1]]
data.tree$tip.label <- c("Neuronal", "Neuronal", "Neuronal","Microglia","Astrocyte",
"Interneuron","Interneuron","Endothelial","Interneuron",
"Oligodendrocytes","Pericyte","Pericyte")
plot.phylo(x = data.tree, font = 1)
i <- c(1,2,3) #Off one bc indexing
nodelabels(pch = 19)
tiplabels(data.tree$tip.label[i], i, adj = 0)

#TSNE coloring just neuronal clusters
experiment.aggregate.0.3@meta.data$neuronal <- "not.neuronal"
wt.cells <- grep("WT", rownames(experiment.aggregate.0.3@meta.data), value = T)
null.cells <- grep("NULL", rownames(experiment.aggregate.0.3@meta.data), value = T)
het.cells <- grep("HET", rownames(experiment.aggregate.0.3@meta.data), value = T)
experiment.aggregate.0.3@meta.data[experiment.aggregate.0.3@meta.data$res.0.3 %in% c(0:2) &
rownames(experiment.aggregate.0.3@meta.data) %in% wt.cells,]$neuronal <- "neuronal.WT"
experiment.aggregate.0.3@meta.data[experiment.aggregate.0.3@meta.data$res.0.3 %in% c(0:2) &
rownames(experiment.aggregate.0.3@meta.data) %in% null.cells,]$neuronal <- "neuronal.NULL"
experiment.aggregate.0.3@meta.data[experiment.aggregate.0.3@meta.data$res.0.3 %in% c(0:2) &
rownames(experiment.aggregate.0.3@meta.data) %in% het.cells,]$neuronal <- "neuronal.HET"
TSNEPlot(object = experiment.aggregate.0.3, group.by="neuronal", pt.size=0.5, do.return = TRUE)

10. Output only neuronal cells
neuronal_cells_to_keep <- rownames(experiment.aggregate.0.3@meta.data[experiment.aggregate.0.3@meta.data$res.0.3 %in% c(0,1,2),])
save(list = c("neuronal_cells_to_keep"), file=paste0(wd, "/", out_dir, "/neuronal_cells.RData"))
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin17.7.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS/LAPACK: /usr/local/Cellar/openblas/0.3.6_1/lib/libopenblasp-r0.3.6.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ape_5.3 scales_1.0.0 cluster_2.1.0
## [4] dynamicTreeCut_1.63-1 gridExtra_2.3 Seurat_2.3.4
## [7] Matrix_1.2-17 cowplot_1.0.0 ggplot2_3.2.1
## [10] knitr_1.25
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.15 colorspace_1.4-1 class_7.3-15
## [4] modeltools_0.2-22 ggridges_0.5.1 mclust_5.4.5
## [7] htmlTable_1.13.2 base64enc_0.1-3 rstudioapi_0.10
## [10] proxy_0.4-23 npsurv_0.4-0 flexmix_2.3-15
## [13] bit64_0.9-7 codetools_0.2-16 splines_3.6.0
## [16] R.methodsS3_1.7.1 lsei_1.2-0 robustbase_0.93-5
## [19] zeallot_0.1.0 Formula_1.2-3 jsonlite_1.6
## [22] ica_1.0-2 kernlab_0.9-27 png_0.1-7
## [25] R.oo_1.22.0 compiler_3.6.0 httr_1.4.1
## [28] backports_1.1.5 assertthat_0.2.1 lazyeval_0.2.2
## [31] lars_1.2 acepack_1.4.1 htmltools_0.4.0
## [34] tools_3.6.0 igraph_1.2.4.1 gtable_0.3.0
## [37] glue_1.3.1 RANN_2.6.1 reshape2_1.4.3
## [40] dplyr_0.8.3 Rcpp_1.0.2 vctrs_0.2.0
## [43] gdata_2.18.0 nlme_3.1-141 iterators_1.0.12
## [46] fpc_2.2-3 gbRd_0.4-11 lmtest_0.9-37
## [49] xfun_0.10 stringr_1.4.0 lifecycle_0.1.0
## [52] irlba_2.3.3 gtools_3.8.1 DEoptimR_1.0-8
## [55] MASS_7.3-51.4 zoo_1.8-6 doSNOW_1.0.18
## [58] parallel_3.6.0 RColorBrewer_1.1-2 yaml_2.2.0
## [61] reticulate_1.13 pbapply_1.4-2 rpart_4.1-15
## [64] segmented_1.0-0 latticeExtra_0.6-28 stringi_1.4.3
## [67] foreach_1.4.7 checkmate_1.9.4 caTools_1.17.1.2
## [70] bibtex_0.4.2 Rdpack_0.11-0 SDMTools_1.1-221.1
## [73] rlang_0.4.0 pkgconfig_2.0.3 dtw_1.21-3
## [76] prabclus_2.3-1 bitops_1.0-6 evaluate_0.14
## [79] lattice_0.20-38 ROCR_1.0-7 purrr_0.3.2
## [82] labeling_0.3 htmlwidgets_1.5.1 bit_1.1-14
## [85] tidyselect_0.2.5 plyr_1.8.4 magrittr_1.5
## [88] R6_2.4.0 snow_0.4-3 gplots_3.0.1.1
## [91] Hmisc_4.2-0 pillar_1.4.2 foreign_0.8-72
## [94] withr_2.1.2 fitdistrplus_1.0-14 mixtools_1.1.0
## [97] survival_2.44-1.1 nnet_7.3-12 tsne_0.1-3
## [100] tibble_2.1.3 crayon_1.3.4 hdf5r_1.3.0
## [103] KernSmooth_2.23-15 rmarkdown_1.16 data.table_1.12.4
## [106] metap_1.1 digest_0.6.21 diptest_0.75-7
## [109] tidyr_1.0.0 R.utils_2.9.0 stats4_3.6.0
## [112] munsell_0.5.0